close all
clear
clc

% 2D FDTD simulation with PLM

%% Common value for simulation
c0 = 299792458; % light speed
e0 = 8.8541878176e-12;% permittivity of free space
u0 = 1.2566370614e-6; % permeability of free space
eta0 = sqrt(u0/e0); % impedance of free space

gigahertz = 1e9;
megahertz = 1e6;

ER_air = 1;
UR_air = 1;
nmax = sqrt(ER_air*UR_air); % refractive index for material, here is air

freq_max = 3.5*gigahertz; % 1GHz

x_length = 2.6; % 1m in x axis
y_length = 2.6; % 1m in y axis


%% Start calculate grid resolution
NRES = 20; % resolve the wave with at least 10 cells, better >=10
lambda_min = c0./freq_max./nmax;
delta_lambda = lambda_min/NRES;

% use free space in this example so material resultion is not needed
%NDRES = 10; % normally 1~4 for resolution of feature size
%delta_d = x_length/NDRES;
%d_temp = min(delta_lambda, delta_d);

d_temp = delta_lambda;

Nx = ceil(x_length/d_temp);
Ny = ceil(y_length/d_temp);
dx = x_length/Nx; % get the size of grid
dy = y_length/Ny; % get the size of grid
d_avg = (dx+dy)/2;

% cells for NPML: NXLO;NXHI;NYLO;NYHI
NPML = [30;30;30;30]; % cells for PML in 1x Grid

%% Start calculate delta_t and tau

Nt = 10; % resuolution for t, normally >=10;
dt1 = 1/Nt/freq_max; 
dt2 = min(dx,dy)/(2*c0); % also consider 1 cell 2 time steps also
dt = min(dt1,dt2);
delta_freq = 1*megahertz; % target frequency resolution, here is 1Mhz

% STEPS>=1/(dt*delta_f) ths is to make sure freq resoltuion in freq domain
% STEPS = 1/dt/delta_freq;

% Pusle width if using Gauss Source
% tau = 1/(2*freq_max);
% t0 = 6*tau;
tau = 3/freq_max; % sin source
t0 = 3*tau;

% another way to cal STEPS is to calculate prop time if freq domain does
% not need to be considered
t_prop = max(nmax*Nx*dx/c0, nmax*Ny*dy/c0);
T_total = 12*tau + 5*t_prop;
STEPS = ceil(T_total/dt); % STEPS calculated by prop time

t_axis=(0:STEPS-1).*dt; % time axis

%% 2x Grid calculation
Nx2 = 2*Nx; Ny2 = 2*Ny;
dx2 = dx/2; dy2 = dy/2;

% calculate ER/UR based on 2x grid
% this example is free space and always 1

ER2 = ones(Nx2,Ny2); % isotropic for free space
UR2 = ones(Nx2,Ny2); % isotropic for free space

% after ER2 and ER2 ready, notmally need to 'smooth' it  by convolution
% this example is free space and no need smooth

% calculate for Ez mode
ERzz = ER2(1:2:Nx2,1:2:Ny2);
URxx = UR2(1:2:Nx2,2:2:Ny2);
URyy = UR2(2:2:Nx2,1:2:Ny2);

%% Gause Source for Ez mode 
Nx_src = [297;302;307;312];
Ny_src = 80;
ER_src = ERzz(Nx_src(1), Ny_src);
UR_src = (URxx(Nx_src(1), Ny_src) + URyy(Nx_src(1), Ny_src))/2;
n_src = sqrt(ER_src(1)*UR_src);

delta_t = n_src*d_avg/(2*c0) + dt/2; % total delay between E and H

A = exp(-((t_axis-t0)/tau).^2);
A(ceil(t0/dt):end) = 1;
Esrc = A.*sin(2*pi*freq_max*t_axis);
% A = -sqrt(ER_src/UR_src); % amplitude of H field
% Esrc = exp(-((t_axis-t0)/tau).^2); % E filed source
% Hsrc = A*exp(-((t_axis-t0+delta_t)/tau).^2); % H field source

%% Calculate PML parameters
sigx = zeros(Nx2,Ny2);
for nx = 1 : 2*NPML(1)
    nx1 = 2*NPML(1) - nx + 1;
    sigx(nx1,:) = (0.5*e0/dt)*(nx/2/NPML(1))^3;
end
for nx = 1 : 2*NPML(2)
    nx1 = Nx2 - 2*NPML(2) + nx;
    sigx(nx1,:) = (0.5*e0/dt)*(nx/2/NPML(2))^3;
end

sigy = zeros(Nx2,Ny2);
for ny = 1 : 2*NPML(3)
    ny1 = 2*NPML(3) - ny + 1;
    sigy(:,ny1) = (0.5*e0/dt)*(ny/2/NPML(3))^3;
end
for ny = 1 : 2*NPML(4)
    ny1 = Ny2 - 2*NPML(4) + ny;
    sigy(:,ny1) = (0.5*e0/dt)*(ny/2/NPML(4))^3;
end

%% fetch PML function onto 1x
% calculate updated coefficients for 2D model
sigHx = sigx(1:2:Nx2,2:2:Ny2);
sigHy = sigy(1:2:Nx2,2:2:Ny2);
mHx0 = (1/dt) + sigHy/(2*e0);
mHx1 = ((1/dt) - sigHy/(2*e0))./mHx0;
mHx2 = - c0./URxx./mHx0;
mHx3 = - (c0*dt/e0) * sigHx./URxx ./ mHx0;

sigHx = sigx(2:2:Nx2,1:2:Ny2);
sigHy = sigy(2:2:Nx2,1:2:Ny2);
mHy0 = (1/dt) + sigHx/(2*e0);
mHy1 = ((1/dt) - sigHx/(2*e0))./mHy0;
mHy2 = - c0./URyy./mHy0;
mHy3 = - (c0*dt/e0) * sigHy./URyy ./ mHy0;

sigDx = sigx(1:2:Nx2,1:2:Ny2);
sigDy = sigy(1:2:Nx2,1:2:Ny2);
mDz0 = (1/dt) + (sigDx + sigDy)/(2*e0) ...
    + sigDx.*sigDy*(dt/4/e0^2);
mDz1 = (1/dt) - (sigDx + sigDy)/(2*e0) ...
    - sigDx.*sigDy*(dt/4/e0^2);
mDz1 = mDz1 ./ mDz0;
mDz2 = c0./mDz0;
mDz4 = - (dt/e0^2)*sigDx.*sigDy./mDz0;

mEz1 = 1./ERzz;

%% Initialize field to zero
Hx = zeros(Nx,Ny); Hy = zeros(Nx,Ny); Ez = zeros(Nx,Ny);
Dz = ERzz.*Ez;
CEx = zeros(Nx,Ny); CEy = zeros(Nx,Ny); CHz = zeros(Nx,Ny);
ICEx = zeros(Nx,Ny); ICEy = zeros(Nx,Ny);
IDz = zeros(Nx,Ny);

%% initialize plot parameters
fig=figure;
set(fig,'Name', 'FDTD 2D PML Simulation');
set(fig,'NumberTitle', 'off');
xa=1:Nx;
ya=1:Ny;
[Y,X] = meshgrid(ya,xa);


%% Core FDTD Loop for 2D Ez mode
for T = 1 : STEPS
    % Calculate curl of E, update CEx
    for nx = 1 : Nx
        for ny = 1 : Ny-1
            CEx(nx,ny) = (Ez(nx,ny+1) - Ez(nx,ny))/dy;
        end
        CEx(nx,Ny) = (0 - Ez(nx,Ny))/dy; % Dirichlet Boundary Conditions
    end
    % Calculate curl of E, update CEy
    for ny = 1 : Ny
        for nx = 1 : Nx-1
            CEy(nx,ny) = -(Ez(nx+1,ny) - Ez(nx,ny))/dx;
        end
        CEy(Nx,ny) = -(0 - Ez(Nx,ny))/dx; % Dirichlet Boundary Conditions
    end
    
    % Update H Integrations
    ICEx = ICEx + CEx;
    ICEy = ICEy + CEy;
    
    % Update H from E
    Hx = mHx1.*Hx + mHx2.*CEx + mHx3.*ICEx;
    Hy = mHy1.*Hy + mHy2.*CEy + mHy3.*ICEy;
    
    
    % Calculate curl of H, update CHz
    CHz(1,1) = (Hy(1,1) - 0)/dx...
        -(Hx(1,1) - 0)/dy;
    for nx = 2:Nx
        CHz(nx,1) = (Hy(nx,1) - Hy(nx-1,1))/dx...
            -(Hx(nx,1) - 0)/dy;
    end
    for ny = 2:Ny
        CHz(1,ny) = (Hy(1,ny) - 0)/dx...
            -(Hx(1,ny) - Hx(1, ny-1))/dy;
        for nx = 2:Nx
            CHz(nx,ny) = (Hy(nx,ny) - Hy(nx-1,ny))/dx...
                -(Hx(nx,ny) - Hx(nx, ny-1))/dy;
        end
    end
    
    % Update D Integrations
    IDz = IDz + Dz;
    
    % Update Dz
    Dz = mDz1.*Dz + mDz2.*CHz + mDz4.*IDz;
    
    % Inject Source
    Dz(Nx_src(1), Ny_src) = Dz(Nx_src(1), Ny_src) + Esrc(T)*ERzz(Nx_src(1),Ny_src);
    Dz(Nx_src(2), Ny_src) = Dz(Nx_src(2), Ny_src) + Esrc(T)*ERzz(Nx_src(2),Ny_src);
    Dz(Nx_src(3), Ny_src) = Dz(Nx_src(3), Ny_src) + Esrc(T)*ERzz(Nx_src(3),Ny_src);
    Dz(Nx_src(4), Ny_src) = Dz(Nx_src(4), Ny_src) + Esrc(T)*ERzz(Nx_src(4),Ny_src);

    % Update Ez
    Ez = mEz1.*Dz;
    
    
    %Plot
    pcolor(X,Y,10*log10(abs(Ez).^2));
    axis ij
    axis equal tight;
    colormap('jet');
    shading interp;
    caxis([-60,0]);
    %colorbar
    hold on
    title(sprintf('Step: %d of %d',T, STEPS));
    
    drawnow;
    hold off
    %pause(0.01);
    
end









